/*
 * matrix.c
 *
 *  Created on: 16/03/2011
 *      Author: Yotam and Gilad
 */
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include "matrix.h"

enum {
	ZERO = 0, INVALID_DIMENSION = -1
};

/*******Internal functions - Not included in the "matrix.h" header file*******/

/**
 * Gets the rows of the matrix.
 * @param matrix The matrix to operate on.
 * @return The rows of the matrix.
 **/
static int getRows(Matrix matrix) {
	return matrix->rows;
}
/**
 * Gets the columns of the matrix.
 * @param matrix The matrix to operate on.
 * @return The columns of the matrix.
 **/
static int getColumns(Matrix matrix) {
	return matrix->columns;
}
/**
 * swap two elements in a matrix.
 * @param matrix The matrix to operate on.
 * @param row The row to change at.
 * @param column The column to change at.
 **/
static void matrixSwapElements(Matrix matrix, int row1, int column1, int row2,
		int column2) {
	Complex temp = matrixGet(matrix, row1, column1);
	matrixSet(matrix, row1, column1, matrixGet(matrix, row2, column2));
	matrixSet(matrix, row2, column2, temp);
}
/**
 * The function checks whether two matrices dimensions are equal.
 * @param matrix1 The first matrix.
 * @param matrix2 The second matrix.
 * @return Will return the result of rows*columns if the dimensions are equal,
 * INVALID_DIMENSION if they are not.
 **/
static int compareDimentions(Matrix matrix1, Matrix matrix2) {
	int rowsMatrix1 = getRows(matrix1), rowsMatrix2 = getRows(matrix2),
			columnsMatrix1 = getColumns(matrix1), columnsMatrix2 = getColumns(
					matrix2);
	if ((rowsMatrix1 != rowsMatrix2) || (columnsMatrix1 != columnsMatrix2)) {
		return INVALID_DIMENSION;
	}
	return (rowsMatrix1 * columnsMatrix1);
}
/**
 * The function gets the result of multiplying between two elements in matrices.
 * @param matrix1 The first matrix.
 * @param row1 The rows of the first matrix.
 * @param column1 The columns of the first matrix.
 * @param matrix2 The second matrix.
 * @param rows2 The rows of the second matrix.
 * @param columns2 The columns of the second matrix.
 * @return Will return the result of multiplying two complex elements.
 **/
static Complex matrixMultiplyElements(Matrix matrix1, int row1, int column1,
		Matrix matrix2, int row2, int column2) {
	return complexMultiply(matrixGet(matrix1, row1, column1), matrixGet(
			matrix2, row2, column2));
}
/**
 * Multiplies one element in a matrix by a scalar.
 * @param matrix The matrix to operate on.
 * @param row The element's row.
 * @param column The element's column.
 * @return The result of the multiplying the element by the scalar.
 **/
static Complex matrixMultiplyElementByScalar(Matrix matrix, int row,
		int column, Complex scalar) {
	return complexMultiply(matrixGet(matrix, row, column), scalar);
}

/**
 * Function gets two matrices, ranks the first matrix and performs the same
 * actions on the second matrix.
 * @param matrix1 - the matrix to be ranked.
 * @param matrix2 - the matrix that we make the action on.
 * @return MATRIX_SUCCESS - if matrix1 became unit-matrix and matrix2 became
 *  inverse matrix .
 * @return MATRIX_SINGULAR - if matrix1 is singular.
 **/
static MatrixResult rankMatrix(Matrix originalMatrix, Matrix matrixToInverse) {
	int columns = getColumns(originalMatrix), rows = getRows(originalMatrix);
	Complex leadingElement = complexCreate(0, 0);
	bool leadingElementFound = false;

	/*for each column of originalMatrix*/
	for (int columnIndex = 0; columnIndex < columns; columnIndex++) {
		leadingElementFound = false;
		/*searching the leading element only in the relevant rows, which are
		 * from the row that equals to the column, down*/
		for (int rowIndex = columnIndex; rowIndex < rows; rowIndex++) {
			leadingElement = matrixGet(originalMatrix, rowIndex, columnIndex);
			if (!complexIsZero(leadingElement)) {
				leadingElementFound = true;
				/*manipulate the row so that the leading element will be
				 * a complex unit*/
				matrixMultiplyRowByScalar(originalMatrix, rowIndex,
						complexInverse(leadingElement));
				matrixMultiplyRowByScalar(matrixToInverse, rowIndex,
						complexInverse(leadingElement));

				/*swap rows so that the leading element will be
				 * located in (i,i)*/
				matrixSwapRows(originalMatrix, rowIndex, columnIndex);
				matrixSwapRows(matrixToInverse, rowIndex, columnIndex);

				/*zero all the other elements in this column*/
				for (int rowToSubtract = 0; rowToSubtract < rows;
						rowToSubtract++) {
					if (rowToSubtract != columnIndex) {
						Complex subtractor = matrixGet(originalMatrix,
								rowToSubtract, columnIndex);
						matrixSubtractRow(originalMatrix, rowToSubtract,
								columnIndex, subtractor);
						matrixSubtractRow(matrixToInverse, rowToSubtract,
								columnIndex, subtractor);
					}
				}
				break;
			}
		}
		/*if a leading element wasn't found in a column*/
		if (leadingElementFound == false) {
			return MATRIX_SINGULAR;
		}
	}
	/*if a leading element wasn't found in a column*/
	if (leadingElementFound == false) {
		return MATRIX_SINGULAR;
	}
	return MATRIX_SUCCESS;
}

/************************End of the internal functions***********************/

Matrix matrixCreate(int rows, int columns) {
	if ((rows <= 0) || (columns <= 0)) {
		return NULL;
	}
	Matrix newMatrix = (struct matrix_t*) malloc(sizeof(struct matrix_t));
	if (newMatrix == NULL) {
		return NULL;
	}
	/**
	 * Allocates space for a complex array with rows*columns elements of type
	 * Complex.
	 **/
	newMatrix->matrixHead = (Complex*) malloc(sizeof(Complex) * rows * columns);
	if (newMatrix->matrixHead == NULL) {
		free(newMatrix);
		return NULL;
	}
	newMatrix->rows = rows;
	newMatrix->columns = columns;
	for (int i = 0; i < (rows * columns); i++) {
		newMatrix->matrixHead[i] = complexCreate(ZERO, ZERO);
	}
	return newMatrix;
}

Matrix matrixCopy(Matrix matrix) {
	if (matrix == NULL) {
		return NULL;
	}
	Matrix newMatrix = matrixCreate(getRows(matrix), getColumns(matrix));
	if (newMatrix == NULL) {
		return NULL;
	}
	for (int i = 0; i < (getRows(matrix) * getColumns(matrix)); i++) {
		newMatrix->matrixHead[i] = matrix->matrixHead[i];
	}
	return newMatrix;
}

Matrix matrixCreateUnit(int size) {
	if (size <= ZERO) {
		return NULL;
	}
	Matrix newMatrix = matrixCreate(size, size);
	if (newMatrix == NULL) {
		return NULL;
	}
	for (int i = 0; i < size; i++) {
		matrixSet(newMatrix, i, i, complexCreate(1, ZERO));
	}
	return newMatrix;
}

void matrixDestroy(Matrix matrix) {
	if (matrix != NULL) {
		free(matrix->matrixHead);
		matrix->matrixHead = NULL;
		free(matrix);
	}
}

Complex matrixGet(Matrix matrix, int row, int column) {
	/** Check if any of the parameters sent to the function is invalid **/
	if ((matrix == NULL) || (matrix->rows <= row) || (row < 0) || (column < 0)
			|| (matrix->columns <= column)) {
		return complexCreateInvalid();
	}
	int numberOfColumns = getColumns(matrix);
	return matrix->matrixHead[column + row * numberOfColumns];
}

MatrixResult matrixSet(Matrix matrix, int row, int column, Complex value) {
	if (matrix == NULL) {
		return MATRIX_NULL_ARGUMENT;
	}
	if ((row < 0) || (matrix->rows <= row) || (column < 0) || (matrix->columns
			<= column)) {
		return MATRIX_OUT_OF_BOUNDS;
	}
	int numberOfColumns = getColumns(matrix);
	matrix->matrixHead[column + row * numberOfColumns] = value;
	return MATRIX_SUCCESS;
}

MatrixResult matrixAdd(Matrix matrix1, Matrix matrix2) {
	if ((matrix1 == NULL) || (matrix2 == NULL)) {
		return MATRIX_NULL_ARGUMENT;
	}
	int matrixDimensions = compareDimentions(matrix1, matrix2);
	if (matrixDimensions == INVALID_DIMENSION) {
		return MATRIX_BAD_DIMENSIONS;
	}
	for (int i = 0; i < matrixDimensions; i++) {
		matrix1->matrixHead[i] = complexAdd(matrix1->matrixHead[i],
				matrix2->matrixHead[i]);
	}
	return MATRIX_SUCCESS;
}

MatrixResult matrixNegate(Matrix matrix) {
	if (matrix == NULL) {
		return MATRIX_NULL_ARGUMENT;
	}
	for (int i = 0; i < ((getRows(matrix)) * (getColumns(matrix))); i++) {
		matrix->matrixHead[i] = complexNegate(matrix->matrixHead[i]);
	}
	return MATRIX_SUCCESS;
}

MatrixResult matrixSubtract(Matrix matrix1, Matrix matrix2) {
	if ((matrix1 == NULL) || (matrix2 == NULL)) {
		return MATRIX_NULL_ARGUMENT;
	}
	int matrixDimensions = compareDimentions(matrix1, matrix2);
	if (matrixDimensions == INVALID_DIMENSION) {
		return MATRIX_BAD_DIMENSIONS;
	}
	for (int i = 0; i < matrixDimensions; i++) {
		matrix1->matrixHead[i] = complexAdd(matrix1->matrixHead[i],
				complexNegate(matrix2->matrixHead[i]));
	}
	return MATRIX_SUCCESS;
}

bool matrixIsEqual(Matrix matrix1, Matrix matrix2) {
	if ((matrix1 == NULL) && (matrix2 == NULL)) {
		return true;
	} else if ((matrix1 == NULL) || (matrix2 == NULL)) {
		//If only one pointer == NULL:
		return false;
	}
	if (compareDimentions(matrix1, matrix2) == INVALID_DIMENSION) {
		return false;
	}
	for (int i = 0; i < getRows(matrix1); i++) {
		for (int j = 0; j < getColumns(matrix1); j++) {
			if (!complexIsEqual(matrixGet(matrix1, i, j), matrixGet(matrix2, i,
					j))) {
				return false;
			}
		}
	}
	return true;
}

MatrixResult matrixMultiply(Matrix matrix1, Matrix matrix2, Matrix* result) {
	if ((matrix1 == NULL) || (matrix2 == NULL) || (result == NULL)) {
		return MATRIX_NULL_ARGUMENT;
	}
	if ((getColumns(matrix1)) != (getRows(matrix2))) {
		return MATRIX_BAD_DIMENSIONS;
	}
	int rows = getRows(matrix1), columns = getColumns(matrix2);
	Matrix newMatrix = matrixCreate(rows, columns);
	if (newMatrix == NULL) {
		return MATRIX_OUT_OF_MEMORY;
	}

	for (int i = 0; i < rows; ++i) {
		for (int j = 0; j < columns; ++j) {
			for (int k = 0; k < getColumns(matrix1); ++k) {
				matrixSet(newMatrix, i, j, complexAdd(
						matrixGet(newMatrix, i, j), matrixMultiplyElements(
								matrix1, i, k, matrix2, k, j)));
			}
		}
	}
	*result = newMatrix;
	return MATRIX_SUCCESS;
}

MatrixResult matrixSwapRows(Matrix matrix, int row1, int row2) {
	if (matrix == NULL) {
		return MATRIX_NULL_ARGUMENT;
	}
	if ((row1 < 0) || (row1 >= getRows(matrix)) || (row2 < 0) || (row2
			>= getRows(matrix))) {
		return MATRIX_OUT_OF_BOUNDS;
	}
	for (int j = 0; j < getColumns(matrix); j++) {
		matrixSwapElements(matrix, row1, j, row2, j);
	}
	return MATRIX_SUCCESS;
}

MatrixResult matrixMultiplyRowByScalar(Matrix matrix, int row, Complex scalar) {
	if (matrix == NULL) {
		return MATRIX_NULL_ARGUMENT;
	}
	if ((row < 0) || (row >= getRows(matrix))) {
		return MATRIX_OUT_OF_BOUNDS;
	}
	for (int j = 0; j < getColumns(matrix); j++) {
		matrixSet(matrix, row, j, (matrixMultiplyElementByScalar(matrix, row,
				j, scalar)));
	}
	return MATRIX_SUCCESS;
}

MatrixResult matrixSubtractRow(Matrix matrix, int row1, int row2,
		Complex scalar) {
	if (matrix == NULL) {
		return MATRIX_NULL_ARGUMENT;
	}
	if ((row1 < 0) || (row1 >= getRows(matrix)) || (row2 < 0) || (row2
			>= getRows(matrix))) {
		return MATRIX_OUT_OF_BOUNDS;
	}
	for (int j = 0; j < getColumns(matrix); j++) {
		matrixSet(matrix, row1, j, (complexSubtract(matrixGet(matrix, row1, j),
				matrixMultiplyElementByScalar(matrix, row2, j, scalar))));
	}
	return MATRIX_SUCCESS;
}

MatrixResult matrixInverse(Matrix matrix, Matrix* result) {
	if ((matrix == NULL) || (result == NULL)) {
		return MATRIX_NULL_ARGUMENT;
	}

	int rows = getRows(matrix);
	int columns = getColumns(matrix);
	if (rows != columns) {
		return MATRIX_BAD_DIMENSIONS;
	}

	Matrix tempMatrix = matrixCopy(matrix);
	Matrix unitMatrix = matrixCreateUnit(rows);
	if ((tempMatrix == NULL) || (unitMatrix == NULL)) {
		matrixDestroy(tempMatrix);
		matrixDestroy(unitMatrix);
		return MATRIX_OUT_OF_MEMORY;
	}
	MatrixResult resultToReturn = rankMatrix(tempMatrix, unitMatrix);
	matrixDestroy(tempMatrix);
	if (resultToReturn == MATRIX_SINGULAR) {
		matrixDestroy(unitMatrix);
		return MATRIX_SINGULAR;
	}
	*result = unitMatrix;
	return MATRIX_SUCCESS;
}

MatrixResult matrixSolveEquationSystem(Matrix A, Matrix b, Matrix* x) {
	if ((A == NULL) || (b == NULL) || (x == NULL)) {
		return MATRIX_NULL_ARGUMENT;
	}

	int rows = getRows(A);
	if ((rows != getColumns(A)) || (rows != getRows(b)) ||
			(getColumns(b) != 1)) {
		return MATRIX_BAD_DIMENSIONS;
	}
	Matrix tempMatrix = matrixCopy(A);
	Matrix solution = matrixCopy(b);
	if ((tempMatrix == NULL) || (solution == NULL)) {
		matrixDestroy(tempMatrix);
		matrixDestroy(solution);
		return MATRIX_OUT_OF_MEMORY;
	}
	MatrixResult resultToReturn = rankMatrix(tempMatrix, solution);

	matrixDestroy(tempMatrix);
	if (resultToReturn == MATRIX_SINGULAR) {
		matrixDestroy(solution);
		return MATRIX_SINGULAR;
	}
	*x = solution;
	return MATRIX_SUCCESS;
}

